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I. INTRODUCTION 



Numerical simulations in fluid dynamics have historically attracted considerable 
attention, through the use of various numerical schemes such as finite-difference and 
finite-element methods. In this study, the direct numerical simulation of a time-dependent, 
axially symmetric, incompressible viscous flow past a rigid sphere was studied by using a 
different numerical scheme based on spectral methods. Spectral methods have become 
increasingly popular in recent years, especially with the development of fast transform 
methods. It has been successfully applied to numerical weather predictions, numerical 
simulation of laminar/turbulent flows and other fundamental problems especially where 
high accuracy is desired for relatively lower memory requirements. 

The drag dependence on Reynolds number was studied and the results were 
compared with the previous results. The stream function vj/ and the vorticity GO were 
visualized through the use of the numerical results and compared with the flow 
visualization pictures for various Reynolds numbers. The boundary-layer thickness on the 
surface of the sphere was also investigated through the use of the results of the 
simulation. The boundary-layer thickness on the surface of the sphere is calculated by 
introducing the power series expansions of the stream function into the boundary-layer 
equations. The results are valid for Re — > analytically. The minimum Reynolds number 

for which the boundary-layer assumption is valid, is investigated by comparing the 
boundary-layer thickness which is obtained from the numerical solution, with analytical 
results in the boundary-layer limit 

It is clear from the foregoing that the objective of this study is to carry out the 
numerical scheme based on the spectral methods, through the use of the vorticity-stream 
function form of the Navier-Stokes equations, on a time-dependent, axisymmetric, viscous 
flow about a rigid sphere. A Matlab program was developed to implement the foregoing 
numerical scheme. The results were compared to the previous studies which were based 
on different methods. However the details of the method itself were not compared to the 
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other methods which were used in previous studies. The expectations are that the results 
will match with the previous results, improve our current ability to compute/predict the 
evolution of flow for a particular geometry and conditions, and, hopefully shed some light 
on the physics of flows which have been heretofore uncalculated. 
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II. BACKGROUND STUDIES 



The collocation method appears to have been first used by Slater (1934) and 
Kantorovic (1934) in specific applications. It was developed as a general method for 
solving ordinary differential equations. In 1938 Lanczos’ work showed for the first time 
that a proper choice of trial functions and distribution of collocation points is crucial to the 
accuracy of solution. This method was revived by Clenshaw (1957), Norton (1963) and 
Wright (1964). The application of Chebyshev polynomial expansions to the initial value 
problems are involved in these studies. 

The spectral collocation method was applied to partial differential equations for 
spatially periodic problems by Kreiss and Oliger (1972) (who called it the Fourier method) 
and Orszag (1972) (who termed it pseudospectral). These were the earliest applications of 
spectral collocation or pseudospectral method. This approach was very attractive because 
of its application to variable-coefficient and even non-linear problems. 

The Galerkin approach which depends on the same trial and test functions, was 
applied to a meteorological model by Silberman (1954). This was the first serious 
application of spectral methods to PDE’s. After Orszag (1969,1970) and Eliasen, 
Machenhauer and Rasmussen (1970) developed transform methods for evaluating 
convolution sums arising from quadratic non-linearity, spectral Galerkin methods only 
became practical for high resolution calculations of such non-linear problems. 

The first unifying mathematical assessment of the theory of spectral methods was 
covered in the monograph by Gottlieb and Orszag (1977). Since then, the theory has been 
extended to cover many different problems, such as variable-coefficient and non-linear 
problems. Applications in fluid dynamics were reviewed in the symposium proceedings 
edited by Voigt, Gottlieb and Hussaini (1984). Other introductory or review articles have 
appeared recently such as Mercier (1981), Hussaini, Salas and Zang (1985), Zang and 
Hussaini (1985), Deville (1984), Gottlieb (1985) and Hussaini and Zang (1987). 
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It is well known that the background of this particular problem of flow over a 
sphere has a long history due to the fundamental interest in external flow over bluff 
bodies. The first theoretical results for steady, uniform flow past a rigid sphere were given 
by Proudman and Pearson (1957) and are valid for Re < 1. For higher Reynolds numbers 
computations have been conducted by Rimon and Cheng (1969), Le Clair, Hamielec and 
Pruppacher (1970) and more recently by Fornberg (1988). All these computations are 
based on finite-difference schemes. Dennis and Walker (1971), Oliver and Chung (1985) 
and Brabston and Keller (1975) used a different computational approach based on an 
expansion of Legendre functions and spherical harmonics. The results of these simulations 
are consistent with the experimental data for the drag forces, and have been used to 
develop improved formulae. Marcus and Tuckerman (1987) developed another technique 
based on the spectral methods and successfully applied it to the simulation of a flow 
between concentric rotating spheres. And more recently Chang and Maxey (1994) have 
used the same method to compute the oscillatory flow about a sphere. 
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III. GOVERNING EQUATIONS 



A. DERIVATION OF THE EQUATIONS 

The well known Navier-Stokes equations of motion of a viscous, incompressible 
fluid are; 



3u _ 1 _ 

— +(u.V)u = — Vp + vV u 
at p 



V.u = 0 



Define dimensionless variables 



* 1 * U T7 * V P T7* 1 T7 

a U F pU 2 a 



where a: radius of the sphere 
U: free stream velocity 
p : density of the fluid 
v : kinematic viscosity of the fluid 
Introduce the dimensionless variables in Eq. (3.1) 

1 U 2 



Uv du U 2 . . 

-+ (u . V)u = 



i 2 8t a 



Vp* + v~y V 2 u* 
pa a 2 



(3.1) 

(3.2) 



(3.3) 



Divide Eq. (3.3) by 



du U.a , „ . U.a . . . 

-r— + (u .V)u Vp +vV u 

3t v v 



(3.4) 



Hereafter the superscript * will be omitted for nondimensional terms. Define the 
Reynolds number based on the diameter, Re d = — — , and take the curl of Eq. (3.4) to 



v 

be able to eliminate the pressure term and introduce vorticity, a> = V x u where the 
vorticity vector, co, is 
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CO — i r (0t+ie(tii9+i<j>C03> 



Because of the absence of any azimuthal fluid motion about the axis of symmetry only the 
^-component of the vorticity survives and the transport equation for the O-component of 
the vorticity is 



|fl) 0 -[Vx(uxfi)^].i f = 



1 



Re d v r sin 0 J 



co. 



(3.5) 



Convert the velocity in Eq. (3.5) by defining the axisymmetric Stream function (\|/) 



u 



= Vx 



¥ 



r sinG 



<|) ; Azimuthal direction 



(3.6) 



In spherical coordinates, Eq. (3.6) is obtained as 



u = 



1 0 ^" 


1 •— 


1 3 \|f" 


_r 2 sin 0 30 . 




_r sin 0 0 r _ 



where 



u = u r i r +u e i e + 0 * 1 * 



(3.7) 



(3.8) 



As noted earlier u,,, =0 in this study. 

In the free stream \ir = ^-Ur 2 sin 2 0 or \ir = -^r 2 sin 2 0 in nondimensional form. 

2 2 

By using the definition of the stream function it can be shown that 

£0j> = 



"a 2 v 1 

* i_ 


1 0/ 


' i avY 


) 


0r 2 r 2 


sin 0 00 \ 


^sin0 00 ) 


r sin 0 j 



(3.9) 



Define an operator D 2 as 



— , 9 2 sin0 9 

D 2 =tt+ 



1 a 



0r 2 r 2 ae 



sin 0 90 . 



(3.10) 
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to obtain the azimuthal component of the vorticity; 

1 =-, 



Ciij. = — 



rsinG 



D > 



(3.11) 



Call [ V x (u x C0n>i<j,)].io in Eq. (3.5) the nonlinear term G(r,9,t) . 

By carrying out the calculations in spherical coordinates it can be shown that 



where 



Vx(uxa>*i*) 



1 



r 3 sin 9 



9(y,D 2 y) 

d(r,p) 



+ 2D 2 \)/L(v)/) 



p = cos9 

p. d 1 d 
L = — r- — + - — 

1-p dr r dp 

d(z,,z 2 ) dz, dz 2 dz, dz 2 

d(x,y) dx dy dy dx 

Combine all terms to get 



(3.12) 



(3.13) 



Re d 1 



l (D2v)+ 2 r 2 



d(V|/,D 2 V| /) _ 



d(r,p) 



+ 2D vyL(\|r) 



= D> 



(3.14) 



Equation (3.14) is called the Vorticity-Stream function form of the Momentum 
equation in axisymmetric spherical coordinate system. For computational purposes, which 
will become clear later, define a modified stream function C(r,9,t), which is related to the 
usual Stokes stream function (\jr) by the relation; 

\y = rCsin9 (3.15) 

Also for any function f (r,9) , define a new operator D 2 such that 

D 2 (fr sin 9) = r sin 9D 2 f 

where 
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D 2 = V 2 



This gives 



where 



1 

r 2 sin 2 0 

g>,=-D 2 C 

u = Vx(Ci,) 



u. = 



1 9 

r sin 0 30 



(Csin0) 



1 9 / _ N 

u ‘ = "7* (rC) 



The modified stream function C is written as the sum C = C + c , 
C = ~ r sin0 , (which is the prescribed uniform-flow free stream condition), 
represents the disturbance produced by the presence of the rigid sphere. 



Returning to Eq. (3.5), rearrange by multiplying both sides with r 



2 

2 ’ 



2 Re a dco 2 Re d 

' = r 

2 9t 2 



r* . ^ =r~ ^^-G(r,0, t) + r 2 D 2 co 



r 2 D 2 C = -r 2 co 



where the nonlinear term 



G(r,0,t) = [Vx(uxt0 fl) i o )].i <t) 



and the operator D 2 is defined as 



D 2 =~t~-| r z — 1 + 



1 



r 2 9r v 9r J r 2 sin 2 0 90 



9 ( A 0 

sin0— 



1 



V 



90 J r 2 sin 2 0 



(3.16) 

(3.17) 

(3.18) 

where 
and c 



(3.19) 

(3.20) 

(3.21) 



(3.22) 
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B. BOUNDARY CONDITIONS 



The flow is determined by numerically solving (3.5) for the vorticity and (3.20) for 
the modified stream function C subject to the boundary conditions; 

u — 0 at r = 1 (3.23) 

u =U(t) “free stream velocity” as r — » (3.24) 

Boundary conditions (3.23) and (3.24) may be reformulated in terms of C as 



dC 

C = 0 , — = 0 on r = 1 
dr 



(0=0, C = — U(t)r sin(0) as r 



(3.25) 

(3.26) 



The Neumann Boundary Conditions in Equation (3.25) are handled with the 
Green’s function method that will be explained later. Homogenous boundary conditions 
on the axis of symmetry (0 = 0, 7t) are automatically satisfied with the choice of the sine 
series expansions in the numerical method as shown later. 

C. BOUNDARY LAYER BEHAVIOR 

The boundary layer thickness on the surface of the sphere can be calculated with 
the method of power series expansion as discussed in, Schlichting (1 987, pp. 235-238). The 
body contour of the sphere with radius a is defined by r(x) = asin(x / a) where x is 
the distance along the surface from the stagnation point. In the boundary-layer limit the 

3 

velocity distribution near the surface of the sphere is given by U(x) = — U„ sin(x / a) . 

2 

By introducing the stream-function into the boundary-layer equations the final 
form of the transformation is obtained for axially symmetrical case; 
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d\\f d 2 \J/ 
dy dxdy 



(3.27) 



fdy 1 dr TT dU d 3 \y 

dx r dx J dy dx dy ' 



with the boundary conditions 
y = 0 



dy 



y = i 



IT = U(x) 

dy 



The potential flow is defined by the series expansion as 



U(x) = u,x + u 3 x 3 +u 5 x s +. 



(3.28a) 

(3.28b) 



(3.29) 



3U„ 1U. 

where u, = ; u, = r- 

1 2 a 3 4 a 3 



T1 



y Sa 
a V v 



and the body contour is expanded to series with the same way; 

r(x) = r 1 x + r 3 x 3 +r 5 x s + (3.30) 

The stream-function is represented by the Blasius series in analogy with (3.29) & 
(3.30) as; 



V(*,y) 



^^{u.xf, (Ti) + 2u 3 x 3 f 3 (ri)+ 



} 



(3.31) 



By introducing equations (3. 29), (3. 30) and (3.31) into equation (3.27) we obtain a 
set of differential equations for the functions f , , f 

If the first two terms of the expanded series are concerned the first equation is 

fr=-f,f,'+-i(f; 2 -l) (3.32) 

where differentiation with respect to r\ is denoted by primes. The boundary conditions are; 

T) = 0 f, = f = 0 (3.33a) 

ri = oo f = 1 (3.33b) 



to 



and the second equation is 




(3.34) 



with boundary conditions 



n = o f 3 = f;=o 



(3.35a) 




(3.35b) 



After solving the Eq. (3.32) and Eq. (3.34) for f, and f 3 the stream function and 

the velocity can be calculated by using Eq. (3.31). In the boundary layer limit the 
boundary-layer thickness (at a given angular location) is chosen to be the radial distance 
from the surface of the sphere where the velocity reaches to its peak value. By using the 
definition of rq it can be shown that 



where r peak is the nondimensional distance where the velocity reaches its peak value. 

By using Eq. (3.36) the calculated values through the use of foregoing 

explanations are compared with the corresponding r peak values which are obtained from 
the numerical results. The Tj^ values are valid as Re d — » «*> and they are compared 
with the boundary-layer thickness obtained from the numerical results. 




(3.36) 
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IV. NUMERICAL REPRESENTATION 



A. COMPUTATIONAL METHOD 

A numerical scheme called the pseudospectral or collocation method based on a 
recent paper by Chang & Maxey(1994) which is in term based on a technique developed 
by Marcus & Tuckerman (1987) is used to solve the Navier-Stokes equations. Spectral 
methods involve representing the solution to a problem as truncated series of known 
functions. From this point of view, spectral methods may be viewed as extreme 
development of the method of weighted residuals (MWR). The trial and the test functions 
are the key elements of the MWR. This is also true for the spectral methods. The key 
difference lies in the choice of the trial and test functions. The trial functions (also called 
the expansion or approximating functions) and the test functions (also known as weight 
functions) are chosen as infinitely differentiable global functions in spectral methods. 
Functions are transformed between physical space and spectral space through the use of 
Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT) which are 
numerically cheap and quick to compute. Although they are not exact due to aliasing and 
truncating errors they rapidly converges and few frequency modes are required to obtain 
accurate results. 

A physical process can be described either in time domain, by the values of same h 
as a function of t, e.g., h(t), or else in frequency domain, where the process is specified by 
giving its amplitude H, as a function of frequency f, that is H(f), with < f < oo. If h(t) is 
real and even then H(f) is real and even or if h(t) is real and odd then H(f) is imaginary and 
odd. The processes can be transformed to either spectral domain or physical domain by 
using the Fourier transform equations. As it will be explained later in this chapter, we are 
dealing with odd functions by using sine series expansions for co, c. The properties of odd 
and even functions are also valid for their derivatives and these properties will help to 
increase the computational efficiency in the calculations. 
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The functions which are introduced into the governing equations are presented 
both in spectral space as a finite series of trial functions and by values at collocation grid 
points in physical space. The vorticity CO and the modified stream function C are expanded 
as a Chebyshev polynomial in the radial direction and as sine series in the azimuthal 
direction. The calculations in the radial direction are done in physical space while the 
derivatives in azimuthal direction are obtained in spectral space. Since each sine term in 
the expansion satisfies the homogenous 0 boundary conditions (to = 0, 'P = 0 at 0 = 0,7t) 
exactly, no further 0 boundary conditions need to be applied. 

In general let a variable, say (0, be expressed as 



co(r,e,t) = Xco n (z,t)Sin6 n 



n=1 



with 



M 

co n (z,t) = £co mn (t)T m (z) 



m=0 



where T m (z) is Chebyshev Polynomial with -1 < z < 1 

n7t 



0 „ = 



" N + l 



z = Cos 



mrc 

M~ 



n = 1,2,.. .N 



m = 0,1, ...M 



(4.1a) 



(4.1b) 



(4.2) 

(4.3) 



An algebraic map is used to map the radial interval 1 < r < u to the interval 
-1 < z < 1 where 






r = 



1 + z V 



1 + L- 

b - z. 



lzl<l 



(4.4) 



with b = 1 H (note that r 0 = 1 , the surface of the sphere) 

r“ - r 0 

L is a scaling parameter which is used to map the radial interval to the z-interval. A 
finite, but large outer r«, was chosen to avoid regularity problems in the radial 
differentiation. Another mapping parameters is defined for higher Reynolds numbers 
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(Re d >20). For small Reynolds numbers a is taken as 1. Two typical grids for two 
different parameters are shown in Figure 1 . 




Figure 1. Typical grids for a=l a)L=8 rinf=250 M=64 N=64 b) L=2 rinf=50 M=64 N=64 



The sine expansions in Eq.(4.1a) and Eq.(4.1b) for (0 and C satisfy the 
homogenous 0 boundary conditions and match exactly the periodicity and symmetry 
conditions in 9. The solution of the elliptic Eq.(3.20) is required to obtain the modified 
stream function C from vorticity © at any time level. Following Marcus & Tuckerman 
(1987) separable derivative operators D r 2 and De 2 are defined as 

r J D 2 =D?+-2-D“ (4.5) 

sm 0 



where 



De 





sin0 — 
30 






sin0 — 
30. 




(4.6) 



In spectral 0-space the effect of operator De 2 at any fixed radial location (physical r-space) 
can be written as 



D 2 .f(0) = 



2fl a 2 . Q _ 3 ' 

sin 0 — — + sm 0 cos 0 1 

30 30 



X f j sin(j0) 






1 



1 



1 



(4.7) 



- (1 + - j 2 ) sin(j0) + - j(l + j) sin(20 + j0) + - j(l - j) sin(20 - j0) 
2 4 4 
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If we write this expression in the form of Dgf(0) = [s ^ Jf j we get 






(i-2)(i-l) 



1 2 ' 



.i = j + 2 



•i = J 



(4-8) 



(i + 2)(i + l) 



,i = j-2 



0 otherwise 



This matrix is NxN if the representation for Dgf is truncated to N Fourier terms. 

Similarly the operation of multiplying a function f (0) by sin" 2 0 produces another NxN 
matrix which is called A matrix where the only non-zero terms are ; 
f al(i) = -i(i + 1) i = j 



A:= = 



i a2(i) = -2i i < j, i + j = even 



(4.9) 



B. TIME INTEGRATION 

Time integration of the Eq.(3.19) and Eq.(3.20) are accomplished through the use 
of an explicit second-order Adams-Bashforth scheme for the nonlinear terms and an 
implicit second-order Crank-Nicholson scheme for the viscous and linear terms. The 
calculations are made in physical r-space and spectral 0-space. If we denote vorticity co(r,t) 
as the vector of N sine series then discretized forms of Eq.(3.19) and Eq.(3.20) are; 



2 Re,] 



co a —co 

t+At t 



At 



= r 



2 Re a ^ an 1 ^ fr^2 ^ A if ® , + “ t+ ., 



[aG, +pG t _ At ] + [D 2 + A] — 



(4.10) 



Rearrange Eq.(4.10) and get ; 



Dj+A-r 2 ^. 

At 



=- 



Dj+A + r 2 ^. 

At 



to, -r ! Re„[aG,-PG,_ 4 ,] (4.11) 



3 1 r 1 

where a = — ; p = - — and G(r,0,t) = |y x (u x co)Ji^ 

jL* £ 



Equation (3.20) is also discretized by using the same method as ; 
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[D r 2 + A]c(r,t + At) = -r 2 co(r,t + At) (4.12) 

The radial derivatives are evaluated by collocation methods and expressed as 
matrix operations on the vector of function values at the grid points in physical r-space. 
For (M+l) collocation points, we introduce the (M+l)x (M+l) matrix operation, that is 

Re, 



D 0) (n) = D 2 + 



a(l) 



At 



' d -r 2 



(4.13) 



where, I is the identity matrix. The operator [D r 2 -r 2 (Re d /At)+A] in Eq.(4.1 1) may now be 



written in block matrix form as 

- D (1) (1) 0 a2(l)I 0 

0 D (1) (2) 0 a2(2)I ... 

0 0 D (I) (3) 0 ... (414) 

0 0 0 D (1) (4) ... 

• ♦ ♦ ••• « « • • • • • • • 

This matrix acts on the vector of a length of (M+l) x N for vorticity coefficients. 
Eq.(4.11) and Eq.(4.12) are upper triangular, block matrix problems in physical r-space 
and spectral 0-space that can be solved by any solver with the Dirichlet boundary 
conditions discussed below. The Eq.(3.32) is treated in the same manner and leads to an 
upper triangular, block matrix problem with each diagonal term on row n replaced by 
D r 2 +al(n)I. 

C. NONLINEAR TERMS 



If one opens the nonlinear term G(r,0, t) in spherical coordinates gets six terms; 



1 3C 
r 3r 30 



1 3o\ 3c 

r 30 3r 



cot(0) 3C 
r “"Sr 



(4.15) 

cot(9) dco^ (o^dC C 3co, 

r 3r r 2 30 r 2 30 

The nonlinear terms in Eq.(4.15) are handled by transforming the spectral 
coefficients to physical space first and transforming the result of the products to spectral 
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space after the multiplication. Basically the derivatives are obtained in spectral space and 
the products are done in physical space and finally the result is transformed to spectral 
space. 

D. GREEN’S FUNCTION METHOD 



In order to enforce the radial boundary conditions in Eq.(4.11) and Eq.(4.12), 
Green’s function method has been developed which gives the correct boundary conditions 
in Eq.(3.25) and Eq.(3.26) and allows the surface vorticity to develop naturally. If CO and 
c consist of two parts, homogenous and particular solutions, we write 

N N 

C0 = (0 P + X^j S j and c=c p + 2>jCj j = 1,2, N (4.16) 

j=i j=i 



If we introduce Eq.(4.16) into Eq.(4.11) and Eq.(4.12), we find the homogenous parts 
with corresponding boundary conditions as follows; 

E 2 63 j (r,0) = 0 (4.17) H 2 ^ =-r 2 63 j (4.18) 

©j(r = 1,0;) = 5jj c j (r = l,0 i ) = O 

c5j(r = °o,0j) = 0 Cj(r = °°,0 i ) = 0 



The particular parts with corresponding boundary conditions are ; 

E 2 O) p (r,0,t + At) = R(r,0, t) (4.19) H 2 c p (r,0,t + At) = -r 2 co p (r,0,t) (4.20) 



co p (r = l,0) = O 
co p (r = ~,0) = O 



c p (r = 1,0) = -Cl r=I 

c (r = <»,0) = O 



E 2 and H 2 are obtained from block matrix defined as in Eq.(4. 14) 

, 3c 9C 

To find A, we enforce — L = — — I . 

dr dr 



If we substitute Eq.(4.16) into Eq.(4.21), we get 



(4.21) 



3r r=1 . 6i 6 3r 



1 X ' = ~T 1 

r=l,e i r=l,e i 



(4.22) 
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Let 



3c dc p 
3r 3r 



Jr=l.e i 



3c. 

= H i and T7 I = A ., 

V* r=I,0j 



(4.23) 



finally from Eq. (4.21 ),Eq. (4.22) and Eq. (4. 23), we find H ; = ^ A t and solve for 



j=i 

A,j as X = A' 1 !!. These X values can now be used in Eq.(4.16) to arrive at the values of go 
and c. 
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V. DISCUSSION OF THE RESULTS 



The results are going to be discussed in three parts. The drag coefficients 
calculated by following the explanation below for various Reynolds numbers are compared 
with the results of Chang & Maxey (1994) and the streamlines contours are compared 
with the experimental flow visualization pictures and finally the boundary-layer behavior in 
the boundary-layer limit is discussed by using the numerical results. 

A. EVALUATION OF THE DRAG 

By integrating the contributions from the stress tensor G- over the sphere surface 

the resultant fluid force on the sphere exerted by the moving fluid is found as; 



where n is the unit outward normal. There is no lift force because the flow is axisymmetric 
and the force F acts parallel to the ffee-stream velocity U. The components of the force 
can be expressed in terms of the spherical polar coordinates as ; 



The no slip boundary conditions in Eq.(3.23) and the condition of incompressible flow in 
Eq.(3.2) causes the dujdr to vanish on the surface of the sphere and the only 
contribution from is through pressure. The shear stress G re is 




(5.1) 




(5.2) 



The normal component of the stress G n is 

9u r 

a rr =-p + 2pv — 



(5-3) 




(5.4) 
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Equation (5.3) and Eq.(5.4) are in dimensional form. By using the non-dimensional 
variables the drag force C d is defined as 

C d =F 1 /7tpa 2 U^ (5.5) 

The drag force C d may be split into two separate parts, C f for the frictional 
component due to shear stress and C p for the pressure component. The frictional 
component is calculated from Eq.(5.2) and Eq.(5.4) as 

C f = --Mo)(r=l,0)sin 2 (0)d0 (5.6) 

Re 0 

The pressure component C p in non-dimensional form 

n 

C p = -J p(r = 1, 0) sin(20)d0 (5.7) 

o 

The pressure variation over the surface of the sphere can be determined once the 
vorticity field is known : 

2 71 9 

p(r = l,0) = p(l,7t)-— f — (rco) I d0' (5.8) 

Ke * or r=i 

A non-dimensional pressure coefficient k(0) can be defined to make the 
calculations easier; 

k(9) = 2p(r = 1,0) (5.9) 

The pressure at the stagnation point, 0 = n , can be evaluated from the radial 
momentum equation as; 

k(rc) = 2p(r = 1 ,jc) = 1 + I dr (5.10) 

Re , r o0 e=Jt 

B. STEADY FLOWS 

In order to verify the numerical method and to provide a comparison, simulations 
for steady, unidirectional flow were carried out for Re <104. The parameters which 
were used in the calculations are listed in Table 1. 
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For all calculations the flow is initially at rest and then a uniform, constant external 
flow is introduced in the direction of 0 = k and the simulation continued till the solution 
has converged to steady state. Two different grids were used in the calculations. A spatial 
discretization of 64x64 grid points was used for Re < 10 and a relatively finer mesh of 
128x128 grid points was used for Re> 20. For Re >20 a mapping parameter “a” is 
defined for finer resolution near the surface. The vorticity diffuses over a large region at 
low Reynolds numbers compared to the sphere radius while it is convected downstream 



with a more complex structure, which needs increased resolution. 



Re 


L 


a 


M 


N 


At 


r _ 
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64 


64 


0.05 
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Table 1. Parameters used in the calculations 



The present results have been summarized in Table 2 with the previous results 
which were calculated by Chang & Maxey (1994). As it can be seen in Table 2 there is a 
general agreement between the present results and the previous results. The separation 
angle, 0 S , and the ratio of bubble length (s) to the diameter of sphere (d) are plotted later 
in this chapter to compare them with the results which were obtained experimentally by 
Taneda (1956). 
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Figure 2 shows the variation of vorticity and streamlines in the vicinity of the 
sphere at Re = 0.1. As it can be noticed the vorticity is diffusing over a large region and 
very slight asymmetry about 9 = 7t / 2 is observable. 





Figure 2. (a) Lines of constant vorticity for Re = 0.1 A© = 0.1 (b) Streamlines A\j/ = 0.3 
Figure 3 shows the constant vorticity and streamlines contours for Re = 5 .The 

asymmetric region is getting more observable about 0 = n / 2 as Reynolds number 

increases. 





Figure 3. (a) Lines of constant vorticity for Re = 5 A© = 0.1 5 (b) Streamlines Ay = 0.25 
In Figure 4 the constant vorticity and streamlines contours are plotted for Re = 40. 
A separation region is clearly observable on both plots. The recirculation in the separation 
region is also noticeable with the asymmetry over the sphere. 
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Figure 4. (a) Lines of constant vorticity for Re = 40 A© = 0.1 (b) Streamlines A\|/ = 0.05 
for ‘ — A\| / = 0.005 for A\J f = 0.00004 for the separation region. 

Figure 5 magnifies the separation region over the constant streamline contours for 
Re = 40. There is a distortion of the vorticity field caused by advection in the flow and a 
separation is clearly visible with a small region of negative streamlines in recirculation. 




Figure 5. Constant streamline contours for Re = 40. A\J/ = 0.00004 in the separation 
region , A\| / = 0.005 for A\}/ = 0.2 for ‘ ‘ 
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Figure 6 shows the constant vorticity and streamlines for Re = 56.5. The 
separation and the recirculation inside the separation region is clearly visible. There is a 
noticeable change in the wake length as the Reynolds number is increased from Re = 40 to 
Re = 56.5. 





Figure 6. (a) Constant vorticity contours for Re = 56.5 ACO = 0.3 (b) Streamline contours 
A\j / =0.005 for ‘ — A\J/ = 0. 15 for A\|/ = 0.0003 for the separation region. 

The constant vorticity and streamline contours for Re = 104 are plotted in Figure 
7. The separation region and recirculation inside the separated flow is clearly visible. The 
increase in the wake length is also noticeable for this particular Reynolds number. It is also 
noticeable from the foregoing plots that the vorticity diffuses to smaller region as 
Reynolds number increases. 





Figure 7. (a) Lines of constant vorticity for Re = 104 Aco = 0.4 (b) Streamlines A\|/ =0.005 
f or * — A\j/= 0.15 for A\J/= 0.0009 for the separation region for Re = 104. 
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Figure 8 magnifies the streamlines for Re = 56.5. The flow visualization picture 
for same Reynolds number is also presented in Figure 9. It is clearly observable that there 
is a similarity between Figure 8 and Figure 9. The wake length, separation location and the 
nature of the streamlines in Figure 8 and Figure 9 look alike in both figures for the 
numerical results and flow visualization picture. 

Figure 10 magnifies the streamlines plot for Re = 104 to present the similarity of it 
with the flow visualization picture in Figure 11 for the same Reynolds number. The 
extension of the recirculating wake to a full diameter in the downstream is clearly visible in 
Figure 10 and Figure 1 1 . The flow is still laminar at this Reynolds number. 

The calculated separation angles and the ratio of the wake lengths to the diameter 
of the sphere are presented in Figure 12 respectively along with the experimental data 
which were obtained by Taneda (1956). There is a good agreement between the present 
results and Taneda’ s data. Taneda has also extrapolated his data and conjectured the 
critical Reynolds number of the “onset of separation” as Re c = 24 .On the other hand, the 

earlier work of Nisi&Porter (1926),which was specially designed to detect the onset of 
separation with ultramicroscopic methods, reported separation at Re = 10 and indicated 
that the onset of separation in an infinite fluid is about Re c = 8 . 

An estimate of Re = 20 is given by the present results as the Reynolds number for 
the onset of separation. This is good agreement with the value of Re c = 20.7 given by 

Chang&Maxey (1994), the value of Re c = 20.5 given by Dennis&Walker (1971) and the 
value of Re c =20given by both Le Clair et al.(1970) and Lin&Lee (1973). Others have 
found separation to first occur at different Reynolds numbers, for example Rimon&Cheng 
(1969) found separation as low as Re c =10. 

Taneda’ s data also suggests that the length of the separation bubble varies linearly 
with the logarithm of the Reynolds number. The present results are consistent with such a 
result and they remark that over this limited range a simple linear relationship with 
Reynolds number is a good approximation. 
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separation region for Re = 56.5 




Figure 9. The flow visualization picture of streamlines for Re = 56.5 ( Archives de 
l’Academie des Sciences de Paris. Payard & Coutanceau 1974 ) 
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Figure 10. Streamlines A\jr =0.005 for ‘ — A\|/ = 0.15 for 0.0009 for the 

separation region for Re = 104. 




Figure 11. Visualization picture for Re = 104 by a thin coating of condensed milk on the 
sphere, which gradually melts and is carried into the stream of water. Taneda 1956. 
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Figure 12.(a) Separation angle 0 S vs. Log Re (b) The ratio of bubble length to diameter 

(s/d) vs. Log Re ‘ Taneda (1956) 



C. BOUNDARY LAYER BEHAVIOR 

Once the stream function was known for the particular flow the velocity over the 
sphere can be calculated by using the Eq.(3.7). The velocity in 0 direction u 9 which is 



parallel to the surface of the sphere can be calculated by using the relation 

1 3\J/ 



u 9 = 



rsin0 3r 

We also know from potential flow theory that the stream function is given by 



(5.11) 



1 



n 






sin 2 0 



(5.12) 



By using the relation in Eq.(5.1 1) u e can be calculated from Eq.(5.12) for potential flow 



as 



u e = 



2r + 1 
2r 3 



sin0 



(5.13) 



Figure 13 shows the u e /sin0 values which are obtained from Eq.(5.13) for 

potential flow along with the numerical results for various Reynolds numbers from this 
study. It is clearly observable that the velocity profile matches with the potential flow at 
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higher Reynolds numbers. The velocity profile for a small Reynolds number is also 
presented in Figure 13 for comparison, and shows that the velocity reaches its free stream 
value over a much larger distance in this Stokes flow regime. 




Figure 13. The u e /sin0 profiles for varying Reynolds numbers (0 = 71 / 4from the nose) 
‘ ‘ represents potential flow, solid lines represent numerical results. 

The approach to boundary-layer behavior is now examined by comparing the 

boundary-layer thickness from numerical results with values of Tjpeak fr° m boundary-layer 

analysis as discussed earlier in Chapter III.C. Figure 14 shows the T| peak values at 



corresponding Reynolds numbers along with the Ti peak value as Re->°o for a particular 
angle. 



The T\ peak values for corresponding Reynolds numbers can be modeled by assuming 
a model such as 



11 = C, 



+ 



_C^ 

Re m 



where C, is the TJ peak value as Re— , C 2 and m are the constants. 



(5.14) 
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By using the T| peak values at Re=30, Re = 56.5 and Re = 104 , C 2 and m can be 
calculated from Eq.(5.14) as C 2 = 2.464 and m = 0.182 to obtain the model as 



T| = C, 



2.464 
f Re 0 182 



(5.15) 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The numerical results have showed a good agreement with the other results which 
were obtained by using different numerical methods. The pressure coefficient, C p , which 

is one of two components of the total drag coefficient, C d , couldn’t be resolved for higher 
Reynolds numbers through the use of present method. But the calculated skin friction 
coefficients, C f , and the stagnation-point pressure k(7t) have been resolved successfully. 

A fine mapping for collocation points is very important to obtain more accurate 
results. The mapping parameters r„ , a and L should be chosen in a way to provide 

enough collocation points in the boundary layer and recirculation regions. A finer mapping 
can increase the accuracy by allowing finer resolution in these large gradient regions. 

The Green’s function method that provides the natural development of vorticity 
along with the Neumann boundary conditions should not be considered as a trivial detail 
and must be used to have a stable convergence process. 

A different FFT technique can be used which does not restrict the number of grid 
points to a power of two. This could reduce the memory requirement by allowing a 
relaxation in the choice of grid points. 

Very good agreement has been obtained with the previous calculations and 
available experimental flow visualization pictures. This study is at a stage where it can 
now be extended to explore higher Reynolds numbers with suitable turbulent models. 
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APPENDIX A. PROGRAM STRUCTURE 



In the first part of the program, problem parameters were defined such that Re, dt, 
L, r_inf, M, N. Calculations which don’t depend on time were also performed in this part. 
In order to gain some additional memory Dl, D2 and D3 matrices in (4.1 1)&(4.12) were 
calculated with the function D_Bmat and saved in terms of diagonal elements only. Non- 
diagonal elements of these matrices were included separately in each iteration. 

Calculations showed that the matrix Ml (derivative operator in r direction) and 
consequently the matrix M3 and the block matrices which are formed by the diagonals of 
Dl, D2 and D3 are ill-conditioned. To avoid the numerical errors in solving the systems of 
equations. Singular Value Decomposition (SVD) method was used. SVD of the matrices 
and calculations of w and c that are used in Green Function method were performed in this 
part. The built-in function SVD in MATLAB, was used for decompositions of the 
matrices in the systems of equations. 

The calculation part basically consists of four steps. In first step, A y matrix was 

calculated from eqn. (4.22) out of the time loop. In second step, inside the time loop the 
particular solutions of w and c were found from eq. (4. 17) and (4. 18). In third step, H that 
is used in Green’s Function method, was calculated and finally homogenous and particular 
solutions are combined to get w and c. 

Most of the calculations were performed in spectral domain. Only nonlinear 
operations were done in physical domain. Nonlinear terms were calculated with the 
functions nonlin for w-c equations . 

A convergence criterion isn’t used in the program to avoid the wrong results which 
are caused by the over-shooting for high Reynolds numbers. Instead, the results are 
checked during the calculations and program was terminated after seeing the steady state. 

Solsvd : This function is used to solve equations after SVD in Matlab. A limit, 
which is defined as a number smaller than the values of S matrix in SVD, is used in this 
function to avoid the numerical errors due to those ill-conditioned matrices. This is 
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avoided by forcing those singular values which are smaller than the defined limit, to 
zero. 10” 14 is used in all calculations. 

Solsvdl : Since this function was used for solving homogenous parts in Green’s 
Function method and performed once before the time-loop, no limit was used in solsvdl. 

Phy : This function was used for transformations of the variables from spectral 
domain to physical domain 

Operat : Derivatives in 0 direction in spectral domain was performed with the 
function, operat. Since these terms were used only in nonlinear terms the output from this 
function is in physical domain. 

L : Mapping parameter 

Re : Reynold’s Number 

dt : Time interval between iterations 

r_inf : Outer boundary 

step : Number of iterations 

N : Number of grid points in 0 direction. 

M : Number of grid points in r direction. 

xcor, ycor : Cartesian coordinates of grid points with respect to center of sphere, 
(used for presentation of the results) 

M 1 : Derivative operator in r direction. 

Dl, D2, D3 : Matrixes defined as in (4.1 1)&(4.12) 
co : Vorticity. 

c : Potential function as in (3.15a)&(3.26) 

wtil, ctil : Homogenous part of the CO and c defined as in Green Function. 

Aij, H, landa : Variables defined as in Green’s function method in chapter IV. 
wp, cp : Particular solution of co and c defined as in Green Function. 

Cf : The frictional component of drag 

kpi : Non-dimensional pressure coefficient at 0 = n 

Cp : Pressure component of drag 
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APPENDIX B. PROGRAM CODES 



%% define program parameters 

% Table 1 list these parameters for every Reynolds numbers 
clear all 

global EMI M2 M3 
L=8 ;Re= 1 ;dt=0 .02; 
r_inf=250;step=10000; 

%% define grid points 



N=32;n=[0:N-l]'; 

teta=n.*pi./N; 

M=32;m=[0:M]'; 
z=cos(pi.*m./M); 
b=l+2*L/(r_inf-l); % map=l 
r=l+L.*(l+z)./(b-z);% map=l 
global N M 
xcor=cos(teta)*r'; 
ycor=sin(teta)*r'; 

E=emat(M); 
disp('E matrix is done') 

%% define Ml, M2, M3 

dz_dr=(b. *(L+r- 1 )-b. *(r- 1 )+L)./(L+r- 1 ). A 2; % map=l 
B=diag(dz_dr); 

M1=B*E; 

R=diag(r);global R Re dt 
M2=R*M1; 

M3=R*(2.*eye(M+l)+M2)*Ml; 

%% define D1 D2 D3 

[D1 ,D2,D3)=D_Bmat(M,N); 

% define initial w c wp cp u 

dum 1=(- 1 .5./r. A 2)*sin(teta’ ))’ ; 

dum2=[dum 1 ;zeros( 1 ,M+ 1 );- 1 .* flipud(dum 1 (2:N,:))]; 

dum3=fft(dum2) ; 
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w=imag(dum3(l:N, :))];w=w’; 
duml=(-0.75+0.25./r. A 2)*sin(teta’))’; 
dum2=[dum 1 ;zeros( 1 ,M+ 1 );- 1 . *flipud(dum 1 (2:N, :))] ; 
dum3=fft(dum2);c=imag(dum3( 1 :N,: ))] ;c=c’ ; 
wp=zeros(M+ 1 ,N);cp=wp; 

% define kmat (for teta der.) and cotteta (for nonlinear terms) 

dum=ones(l,M+l); 

kmat=n*dum; 

cotteta=[0;cot(teta(2:N))]*dum; 
global kmat cotteta teta 

% svd decomposition for D1 D3 and define rsquare (nonlinear) 

Dlu=zeros(N*(M+l),M+l); 

D 1 1=D 1 u;D 1 p=D 1 u;D31=D 1 u ;D3u=D 1 u;D3p=D 1 u; 
rsquare=zeros(M+ 1 ,N); 
for i=l:N 

bas=(i- 1 )*(M+ 1 )+ 1 ; 

son=bas+M; 

rsquare(:,i)=r. A 2; 

Dl(bas,:)=[l zeros(l,M)]; % modify D1 for w @ r=inf 
Dl(son,:)=[zeros(l,M) 1]; % modify D1 for w @ r=l 
[Dll(bas:son,:),Dlu(bas:son,:),Dlp(bas:son,:)]=svd(Dl(bas:son,:)); 
D3(bas,:)=[l zeros(l,M)]; % modify D3 for c @ r=inf 
D3(son,:)=[zeros(l,M) 1]; % modify D3 for c @ r=l 

[D31(bas:son,:),D3u(bas:son,:),D3p(bas:son,:)]=svd(D3(bas:son,:))', 

end 

global rsquare 

% define c_ c_teta c_r u_ u_teta u_r 

c_=(0.5.*sin(teta))*r’; 
c_teta=0.5.*cos(teta))*r'; 
c_r=-(0.5*sin(teta))*ones( 1 ,M+ 1 ); 
global c_ c_teta c_r 

% define c boundary at r=l 

dum4=-0.5 . *si n(teta) ; 

dum5=[dum4;0;- 1 .* flipud(dum4(2:N,:))]; 

dum6=fft(dum5); 

cbr 1 =imag(dum6( 1 ,N))’ ; 
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% begin green function wtil ctil 

wtil=zeros(M+ 1 ,N A 2) ; 
ctil=wtil; 



forjj=l:N 
RR=zeros(M+ 1 ,N); 
duml=[zeros(l,jj-l) 1 zeros(l,N-jj)]; 
dum2=[duml 0 -l.*fliplr(duml(2:N))]; 
dum3=fft(dum2); 

dum4=imag(dum3(l:N));if jj==l;dum4=real(dum3(l:N));end; 
RR(M+l,:)=dum4; 

wtil(:,jj*N)=solsvdl(Dll((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),. 

Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:)*— 

RR(:,N)); 



for i=N-l:-l:l 
sum=zeros(M+ 1,1); 
ii=(jj-l)*N+i; 
for j=i+l:N 

if rem(i+j,2)=2 

sum=sum+(-2*(i-l)).*wtil(:,(jj-l)*N+j); 

end 

sum( 1 )=0; 

sum(M+l)=0; 

end 

RR_=RR(:,i)-sum; 

wtil(:,ii)=solsvdl(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M, 

Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),„. 

Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),„. 

RR_); 



tempwtil=(- 1 .*rsquare).* wtil(:,(jj- 1)*N+1 :jj*N); 
tempwtil(M+ 1 , :)=zeros( 1 ,N) ; 
tempwtil( 1 ,:)=zeros( 1 ,N); 

ctil(:,jj*N)=solsvdl(D31((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,;),.. 

D3u((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

D3p((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

tempwtil(:,N)); 
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for i=N- 1 1 : 1 
ii=(jj-l)*N+i; 
sum=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)=2 

sum=sum+(-2*(i-l)).*ctil(:,(jj-l)*N+j); 

end 

sum(l)=0; 

sum(M+l)=0; 

end 

wtil_=tempwtil(;,i)-sum; 

ctil(:,ii)=solsvd 1 (D31((i- 1 )*(M+ 1 )+ 1 :(i- 1 )*(M+ 1 )+ 1 +M, 
D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M, 
D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
wtilj; 
end 
end 

% find Aij matrix 



for jj=l:N 

blok=ctil(: ,(jj- 1 ) *N+ 1 :jj *N); 
duml=phy(blok')'; 
dum2=(M 1 *dum 1 )’; 

dum3=[dum2; zeros(l,M+l) ; -l.*flipud(dum2(2:N,:))]; 
dum4=fft(dum3);dum5=imag(dum4(l:N,:));dum5=dum5'; 
ctilr_ 1 (jj, 1 :N)=dum5(M+ 1 
end 

Aij=ctilr_l'; 

% find svd of Aij for landa solving 

Aij(l,:)=[l zeros(l,N-l)]; 

[Aiju,Aijs,Aijv]=svd(Aij); 



% begin time loop 
for ops=l:step 

if ops=l ;prc=c;prw=w;pru=u;end 
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% find R1 



Rl=zeros(M+l,N); 

Rl(:,N)=D2((Nd)*(M+l)4d:(N-l)*(M+l)+l+M,:)*w(:,N); 
for i=N-l:-l:l 
summ=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)=2 
summ=summ+(-2*(i- 1)).* w(:,j); 
end 

end 

R 1 ( : ,i)=D2((i- 1 )*(M+ 1 )+ 1 : (i- 1 ) * (M+ 1 )+ 1 +M, :)*■ w( : ,i)+summ; 
end 

R1=-1.*R1; 

% find R2 

R2a=nonlin(w,c); 

R2b=nonlin(prw,prc); 

R2=3/2.*R2a- l/2.*R2b; 

R2=(- 1 *Re).*rsquare.*R2; 

Rtot=Rl+R2; 

prw=w;prc=c; 



%% green function solution 
% find new wp 
%% apply BC for wp 

Rtot( 1 ,:)=zeros( 1 ,N); % BC for wp=0 at r_inf 
Rtot(M+l ,:)=zeros(l ,N); % BC for wp=0 at r_l 

wp(: ,N)=solsvd(D 1 1((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 +M,:) t ... 
Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

D 1 p((N- 1 ) *(M+ 1 )+ 1 :(N- 1 )* (M+ 1 )+ 1 +M,:), ... 
Rtot(:,N)); 



for i=N-l:-l:l 
sum=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)=2 
sum=sum+(-2*(i-l)). :,: wp(:,j); 
end 
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sum(l)=0; 

sum(M+l)=0; 

end 

Rtot_=Rtot(:,i)-sum; 

wp(:,i)=solsvd(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M, 

Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

Rtot_); 



end 

% solve for new cp 
tempwp=- 1 .*rsquare.*wp; 

tempwp(M+l,:)=cbrl; % BC for cp=-c_r at r_l zero for pure spin 
tempwp(l,:)=zeros(l,N); % BC for cp=0 at r_inf 

cp(:,N)=solsvd(D31((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),.„ 

D3u((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

D3p((N-l)*(M+l)+l:(N-l)*(M+l)+l+M, 

tempwp(:,N)); 



for i=N-l:-l:l 
sum=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)=2 
sum=sum+(-2*(i-l)).*cp(:,j); 
end 

sum(l)=0; 

sum(M+l)=0; 

end 

w_=tempwp(: ,i)-sum; 

cp(:,i)=solsvd(D31((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M, 

D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M, 

w_); 



% find landa 
% find H 

duml=phy(cp')’; 
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dum2=(Ml*duml)'; % NxM+1 
dum3=- 1 . *c_r-dum2 ; 

dum4=[dum3 ;zeros( 1 ,M+ 1 ) ;- 1 . *flipud(dum3(2 :N,:))] ; 

dum5=fft(dum4); 

dum6=imag(dum5(l :N,:)); 

dum6=dum6'; % M+lxN 

H=dum6(M+ 1 ,:);H=H'; % Nxl 
H(1)=0; 

landa=solsvd(Aiju,Aijs,Aijv,H); 

% find complete w & c 

dum 1 c=zeros(M+ 1 ,N) ;dum 1 w=dum 1 c ; 
for jj=l:N 

blokc^til(:,(jj-l)*N+l:jj*N); 
blokw=wtil(:,(jj-l)*N+l:jj*N); 
dum2c=landa(jj).*blokc; 
dum2w=landa(jj).*blokw; 
dum 1 c=dum 1 c+dum2c ; 
dum 1 w=dum 1 w+dum2w; 
end 

w=wp+dumlw; 

c=cp+dumlc; 

res=phy(w'); 

resp=phy(prw’); 

resc=phy(c’); 

resc=resc+c_; 

% Calculate Cf for every iteration 
dum 1 =res(: ,M+ 1 ); 
dum2=(sin(teta)). A 2; 
dum3=dum 1 .*dum2; 
cf(ops)=-4/Re*trapz(teta,dum3); 

% Calculate k(7t) for every iteration 

duml=w’; 
for i=l:M+l 
sum(i)=0; 
for j=l:N 

f=((j- 1 )*duml(i,j).*(- 1 . A (j- 1)))-/N ; 
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sum(i)=sum(i)+f; 

end 

end 

dum2=( 1 ./r).*sum; 

kpi(ops)=( 1 +(8/Re)*trapz(r,dum2)); 

% Check for convergence 

if max(max(isnan(w)))==l ;error(' it blowed up! ! ! ! ');end 
% Save the variables to a file 
save tez.mat c w Re ops r_inf cf kpi L N M dt; 
end % end of time loop 

%%%%% 

function [Dl,D2,D3]=D_Bmat(M,N) 

global M3 R dt Re 

k=(M+l)*N; 

Dl=zeros(k,M+l); 

D2=zeros(k,M+ 1 ) ; 

D3=zeros(k,M+ 1 ) ; 
for i=l:N 
for j=i:N 
if i=j 

Dl((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 

=M3+(- 1 *(i- 1 )*(i). *eye(M+ 1 )-Re/dt *R A 2); 

D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 

=M3+(- 1 *(i- 1 )*(i).*eye(M+l)+Re/dt.*R. A 2); 

=M3+(- 1 *(i- 1 )*(i) *eye(M+l )); 
end 
end 
end 

%%%%%%% 

function E=emat(M) 

E=zeros(M,M); 

m=0:M; 

z=cos(pi.*m./M); 
for i=0:M 
for j=0:M 
if i~=j 

if ((i==0 I i==M) & (j==0 I j==M)) I ((i~=0 & i~=M) & (j~=0 & j~=M)) 
c=l;else 

if (i==0 I i==M) & (j~=0 & j~=M) 
c=2;else 
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if (i~=0 & i~=M) & (j==0 1 j==M) 
c=l/2;end;end;end; 

E(i+ 1 ,j+ 1 )=c*(- 1 ) A (i+j)/(z(i+ 1 )-z(j+ 1 )); 
else; 
if i==0 

E(i+1 ,j+l )=(2*M A 2+1 )/6; 
else 
if i=M 

E(i+ 1 j+ 1 )=(2*M A 2+l)/(-6); 
else 

E(i+ 1 ,j+ 1 )=- 1 *z(j+ 1 )/(2*( 1 -z(j+ 1 ) A 2)); 
end 
end 
end 
end 
end 

%%%%%%% 
function R2=nonlin(w,c) 

global kmat rsquare cotteta N M Ml c_ c_teta c_r teta 
% variable list 

% 

% dw/dteta=nonl dc/dteta=non2 dw/dr=non3 dc/dr=non4 (in phy) 

% 

% transpose matrices for column operations in teta direction 

r=sqrt(rsquare)'; 

w=w';c=c'; 

nonl=operat(w); 

non2=operat(c); 

wph=phy(w);cph=phy(c); 

non3=[Ml*wph']'; 

non4=[Ml*cph’]';non4(:,M+l)=-0.5.*sin(teta); 

term 1 =- 1 Jr. *non3.* (non2+c_teta); 

term2= 1 ./r. *non 1 . * (non4+c_r) ; 

term3=-l.*cotteta./r.*wph.*(non4+c_r); 

term4=- 1 . * cotteta./ r. * non3 . * (c_+cph) ; 

term5=wph./ rsquare’ . * (non2+c_teta) ; 

term6=(c_+cph)./rsquare'.*non 1 ; 

termO=term 1 +term2+term3+term4+term5 +term6 ; 

ter=[termO;zeros(l,M+l);-l.*flipud(termO(2:N,:))]; 

RR2=fft(ter); 

R2=imag(RR2( 1 :N,:));R2=R2'; 
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%%%%%%%%% 



function y=operat(x) 
global N M kmat 
re=zeros(N,M+ 1 ); 

X=re+i.*x; 

Y=i.*kmat.*X; 

Y_ful=[Y ;zeros( 1 ,M+ 1 );flipud(Y(2:N,:))] ; 

y=ifft(Y_ful); 

y=real(y(l:N,:)); 

% % % % % % % % % 

function y=phy(x) 
global N M 
re=zeros(size(x)); 

X=re+i.*x; 

Y_ful=[X;zeros( 1 ,M+l);conj(flipud(X(2:N, :)))]; 

dum=ifft(Y_ful); 

y=real(dum(l :N,:)); 



%%%%%%% 



function y=solsvdl(u,s,v,b) 
w=diag(s); 

wmin=max(w)* 1 e- 1 2; 
for i=l:length(w) 

if w(i)<wmin;ww(i)=0;else;ww(i)=l/w(i);end 
end 

www=diag(ww); 

y=v*www*(u'*b); 



%%%%%%% 

function y=solsvd(u,s,v,b) 
w=diag(s); 

%wmin=max(w)*le-8; 
for i=l:length(w) 

% if w(i)<wmin;ww(i)=0;else;ww(i)=l/w(i);end 
ww(i)=l/w(i); 
end 

www=diag(ww); 

y=v*www*(u'*b); 
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